Quantum Liouville-space trajectories for dissipative systems 



O 

o 



> 

in 
in 
o 

O 



^ • 
Oh! 

a ■ 



X 



Jeremy B. MaddoxQ and Eric R. Bittneij] 

Department of Chemistry, University of Houston, Houston Texas 77204 
(Dated: February 1, 2008) 

In this paper we present a new quantum-trajectory based treatment of quantum dynamics suitable 
for dissipative systems. Starting from a de Broglie/Bohm-like representation of the quantum density 
matrix, we derive and define quantum equations-of-motion for Liouville-space trajectories for a 
generalized system coupled to a dissipative environment. Our theory includes a vector potential 
which mixes forward and backwards propagating components and non-local quantum potential which 
continuously produces coherences in the system. These trajectories are then used to propagate an 
adaptive Lagrangian grid which carries the density matrix, p{x,y), and the action, A{x,y), thereby 
providing a complete hydrodynamic-like description of the dynamics. 

PACS numbers: 03.65.x 05.30 34.10-fx 



The causal or hydrodynamic trajectory model of quan- 
tum mechanics provides a useful glimpse into the dy- 
namics underlying the quantum wavefunction.|Q Start- 
ing from the time-dependent Schrodinger equation, one 
can derive quantum trajectory equations by casting the 
wavefunction in polar form, ip — Rexp(iS/h) and sub- 
sequently separating ip into real and imaginary compo- 
nents. The resulting equations are a continuity equa- 
tion for the quantum density, p = , and "Newtonian" 
equations for trajectories. The distinguishing feature of 
these quantum equations-of-motion is the presence of 
the non-local quantum potential, Q, which both intro- 
duces non-local coupling between the trajectories and 
represents a differential geometric constraint between the 
extrinsic curvature invariants of a surface generated by 
z = C\nR and the action per unit volume of a trajectory 
element. This prescription introduced independently by 
de Brogliej^ and Madelung[^, and later developed by 
BohmM has been used primarily as an interpretive tool 
where one first solves the time-dependent Schrodinger 
equation for the wavefunction, and then uses ip 

as an ancillary field to drive an ensemble of quantum 
trajectories. From this construction one can develop in- 
terpretive models of barrier tunneling, arrival times, and 
various other quantum effects based upon these trajecto- 
ries. Furthermore, there has been a recent surge of ac- 
tivity to use this approach to develop new adaptive-grid 
based methods for performing quantum dynamical cal- 
culations. Such ab initio approaches construct ip directly 
from the dynamical information obtained in computing 
the trajectories. [|[ ||, @, |§1 While it is unclear at this 
time whether such approaches offer considerable compu- 
tational advantages over more standard finite-basis set or 
fixed grid approaches, the appeal is that by moving to a 
particle-based description, computational effort required 
scales almost linearly with the number of trajectories in- 
troduced. 



* email :Jmaddox@bavou. uh.ed u 
t email :|bittncr@uh.cdu 



In this letter we derive a causial-trajectory approach 
for the quantum density matrix under dissipative condi- 
tions. So far as we know, this is the first presentation 
of such an approach and certainly is the first ab initio 
application. The formal construction presented here is 
similar in spirit to the Wigner phase-space representation 
which is widely used for quantum dissipative dynamics 
in both the quantum and semi-classical regimes; however, 
our approach does not require the definition of a quan- 
tum phase space and retains all of the physical properties 
associated with the quantum density matrix. Secondly, 
we describe an ab initio approach for computing the tra- 
jectories based upon a finite-element/least squares pro- 
ceedure. This provides a highly efficient and adaptive 
scheme which can be applied to a wide variety of sys- 
tems. Moreover, the computational effort scales almost 
linearly with number of finite-elements/trajectories. In 
this letter we focus upon the construction of the under- 
lying theory with application towards the relaxation of a 
tagged oscillator in contact with an environment. 

The theory is initiated by writing the density ma- 
trix as p{x,y) — ip{x)'ip* (y)- If p represents the en- 
tire system-l-bath ensemble, its evolution is given by the 
Liouville-von Neumann equation, p = Lp, where L is the 
full Liouville operator which we decompose into system 
and bath components: L = Ls + Ljj + Lsb- If we take the 
bath as an ensemble of harmonic oscillators with masses, 
TO„, and frequencies, LOm which is at thermal equilibrium 
at t = coupled to the system coordinate, x, via a linear 
combination x — c„g„, then the response of the bath 
to this coupling is encoded in the correlation functions 
fi{v{t) + irj{t)) ~ {x{t)x{0)) / T where (. . .)t denotes the 
expectation value taken in the thermal equilibrium state 
of the bath. The real and imaginary parts of this correla- 
tion function are the noise and dissipation kernels. This 
procedure is starting point of many quantum statistical 
mechanical treatments and the result is either a master 
equation for ps, or a path- integral treatment for propaga- 
tor. Various forms of the master equation have been pre- 
sented in the literature by a number of groups including 
Caldeira and Legg ett jlH , Unruh and ZurekQ, and by 
Hu and co- workers Each of these is slightly different 
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in final form due to the various approximations made over 
the course of the derivation. Consequently, each posses 
certain regions of validity. Without loss of generality, we 
will focus our attention here on the Caldeira-Leggett [p] 
(CL) model whereby the noise and dissipation terms take 
the form i^{t) = -/h/X^Sit) and r]{t) = --fmS'{t). This 
model is strictly valid in the high-temperature regime 
where the thermal excitations in the bath dominate and 
fluctuations of the vacuum can be ignored. Here, the 
master equation for the reduced density matrix is given 
by 



dtp = Lsp - j{x - y){dxP - dyp) 
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{x-yfp. (1) 



Here, Lg denotes the system Liouvillian where the po- 
tential energy term has been renormalized due to the 
interaction with the bath. For a harmonic system, the 
renormalized oscillator frequency is related to the bare 



frequency, oj, by 
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2juj/t:. The last two terms 



are related to interaction between the system and the 
bath. The second being responsible for population re- 
laxation and the third for coherence relaxation. Finally 
A = h/y^2mkT is the thermal de Broglie wavelength and 
7 is the coupling strength. 

To make the connection to a quantum trajectory 
scheme, we write the density matrix in complex polar 
form, ps{x,y) = exp{g(xTy) + iA{x,y) /fi), and insert 
this into the master equation for ps{x,y). In what fol- 
lows, we shall take x as the coordinate associated with 
ip{x) and y as the coordinate associated with the complex 
conjugate, ij}*{y). Our Liouville space is then the two- 
dimensional configuration space {x, y} with phase-space 
defined as {x,px, y^Py} where Px — dxA and Py = —dyA. 
The — sign in the definition of the canonical momentum 
in y reflects the time-reversed dynamics in the y direc- 
tion. In a geometric-optical construction, the wavefront 
for the density matrix is given by contours of constant 
S{x) + S{y) = c with velocities and momenta are both 
normal to these curves. For non-dissipative dynamics, 
the Hamilton- Jacobi equation for the action A{x,y) is 
separable into forward (in x) and backward (in y) prop- 
agating components: 

OA 1 , , 



dt 

where V — ^(a;) - 
quantum potential. 
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[Px - Pv 



V-Q, 



V{y) and Q - Q{x) - 
Qix) = -h\dlg + dxi 



(2) 

Q{y). The 
dx9)/2m, in- 



troduces non-local interactions between the trajectories 
but does not couple motion in x and y. The resulting 
equations of motion take the form, mv^ = p^, where 
upon taking the time derivative yields Bohm's quantum 
equations of motion in x and y: mi^ — ±dp,{V + Q) {— 
for X and -I- for y). We also have the continuity equation 
for g: 2g ~ —d^v^. 

Taking the CL equation and using the substitution 
logp = g + iA/h a,s before again yields a Hamilton- Jacobi 
equation for A(x,y): 



dtA 
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{pi - pI) - l{x ~y){px-py) + V+Q^ 0. 




FIG. 1: Diagonal trajectories for underdamped (A) and crit- 
ically damped (B) cases at kT — 2.5hui. Dashed line corre- 
sponds to a semi-classical approximation for (^(t)) (c.f.Eq. ^) 



However, this equation now sports an explicit coupling 
between the forward and backward propagating paths. 
This dissipative coupling term involving j(x — y) can be 
cast as a vector potential, F^ = ±7(x— ?/) where — is for x 
and -|- is for y. Using this we define the components of the 
material velocities as mi^ = P/i + T/i which describe the 
velocities of the particles in the presence of dissipation. 

At this point we find it easier to work in the rotated 
coordinate frame given by ^ = (x -I- y)/-\/2 and t] = {y — 
/2. In this frame, the quantum Lagrangian is 



C = —mS^fj -\- 2771777^ — Q — V. 
Here the quantum potential is given by 

Q = — (5^a,,5 + %a,,5r), 



(3) 



(4) 



and V — V{x) — V{ii). In this new frame, the canonical 
momenta are given by = —dC/dfj and p,, = —dC/dS^. 
This modified relation between the canonical momenta 
and the Lagrangian is again due to the time-reversed dy- 
namics in y. The resulting material velocities are then 



77 — Pri/m + 2^7] and ^ ~ p^jm. 



(5) 



Finally, the continuity equation for g is given in the ro- 
tated frame as 



dt 



+ 7- 



27 2 



(6) 
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This equation includes a term proportional to rf' which 
dampens components of the density matrix which are far 
from the diagonal axis leaving the elements along the ^- 
axis unchanged. This is the process of decoherence, i.e. 
the dynamical diagonalization of the density matrix due 
to the interaction with the environment. |14[] 

We take as an example problem the case of harmonic 
oscillator embedded in a dissipative bath. Such a sys- 
tem could easily represent a vibrational coordinate of a 
molecule that has been displaced from its equilibrium 
position by some action at < = 0. Our computational 
methodology, is based upon the moving weighted least- 
squares (MWLS) method reported by Lopreore and Wy- 
att and reader is referred to Ref . for details of this ap- 
proach and to Ref. for our implementation as applied 
to quantum wavepacket dynamics. | [T5[ | 

Taking k as the renornialized force constant and lo^ ~ 
k/m^ the quantum Lagrangian is 

C = —m^i) + 2m7^?7 + k^r/ — Q. (7) 

Using the Euler-Lagrange equations, we can derive equa- 
tions of motion for the Liouville space quantum trajec- 
tories: 

^ = -c^'e --Qn- 27C (8) 

TO 

= j_Q , 2777. (9) 

TO 

Defining the semi-classical limit as ?i ^ so that Q van- 
ishes, the equations of motion are trivial to solve. For 
motion in ^ we have damped harmonic motion which is 
easy to understand in terms of pure populational relax- 
ation. However, in i], the trajectories diverge towards 
±00. 

This introduces a small problem in our computational 
scheme, namely, all of the non-diagonal grid points are 
swept away as the calculation progresses. To circumvent 
this, we re-meshed the grid periodically by discarding the 
off-diagonal trajectories and re-initiating p, g, and A at 
the new points via least-squares interpolation. We also 
took advantage of the parity of these functions under 
77 — > —rj by explicitly propagating only the upper trian- 
gle of each and selecting polynomial bases which carried 
the appropriate symmetry. To whit: g{S,,'i]) = gi£,,—v)j 
Ai^,rj) = -A{t-7j), Qitv) = -Qii,-V), vdtv) = 
'^(.{^1 ^^r;(Ci '7) = ~''^ri(^^ ^v)- Thcsc Symmetries were 
enforced throughout the calculation. 

As an initial example, we consider the relaxation of 
a density matrix corresponding to the ground state dis- 
placed from equilibrium by a shift in x. 

9iO)^9o-:^i{^-ior+rf) (10) 

In the plots shown here, we scaled the time units in terms 
of the oscillator period r = 2tt/uj and lengths by the 
initial width of the gaussian density matrix, = h/muj. 

In the Fig. 1, we show an ensemble of quantum trajec- 
tories corresponding to diagonal element of the density 



matrix (i.e. ij = 0) at kT = 2.5?iu! in the underdamped 
{j/uj = 0.25) and critically damped (7/w = 1) regimes. 
The superimposed dashed line is the semi-classical tra- 
jectory obtained by solving Eq.||in the limit of ?i ^ 0. As 
one expects, the the semi-classical result follows the peak 
of the gaussian. Furthermore, the diagonal trajectories 
are not allowed to cross one another. This is a directly 
analogous to the non-crossing rule for ordinary quantum 
trajectories. Each diagonal element carries a unique pop- 
ulation trajectory which must remain single-valued. As 
time progresses the trajectories become stationary corre- 
sponding to thermal equilibrium. 

In Fig. 2 we show the time evolution of the gaus- 
sian coefficients cr^ and cr,, for various temperatures and 
coupling constants. At thermal equilibrium, these quan- 
tities should become stationary and can be computed 
from equilibrium statistical mechanics, 

4.^9 = (7lcoth{hujP/2) (11) 
al^^ = altanh{hcu/3/2). (12) 

For the critically damped case, we see exponential re- 
laxation back to equilibrium. For ctj, the relaxation is 
largely independent of temperature and each case reaches 
the equilibrium value at t/r = 1. The relaxation of the 
coherence width, cr^, shows a strong temperature depen- 
dence, relaxing faster at higher temperatures. Moreover, 
the time-scale for coherence relaxation is considerably 
shorter than in cr^. This is due to the fact that the ther- 
mal de Broglie wavelength. A, sets both a length and 
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FIG. 2: Relaxation of cr^ and a,, at various temperatures in 
the overdamped (7/0; = 1.75), underdamped (7/cj = 0.25), 
and critically damped (7/0; = 1) regimes. 
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time-scale for decoherence. pq As temperature increases, 
this length decreases and the relaxation rate increases. 
The underdamped case shows rather interesting dynam- 
ics by relaxing through a series of intermediate plateaus. 
Finally, the overdamped case relaxes slower than the crit- 
ically damped without the plateaus seen in the under- 
damped case. 

The picture we are lead to is that the dissipative cou- 
pling to the environment causes a net flux of trajectory el- 
ements (representing population coherence information) 
towards rj — *■ ±(X) and causes the populations along the 
r] — axis to relax to some lowest energy configuration. 
At T = this would be the quantum ground state of 
the system. At finite temperature the system becomes 
stretched in ^ reflecting a thermal distribution of energy 
states. Furthermore, the coherence length as set by the 
de Broglie wavelength, becomes more and more short 
ranged as T increases causing the system to become local- 
ized in 77, effectively diagonalizing the density matrix. In 
the equations-of-motion, both the quantum potential and 
the vector potential, F, force particles to stream outwards 
in r] away from the diagonal. Consequently, even when 



the system becomes stationary (i.e. g = 0) the trajecto- 
ries themselves remain in constant motion reflecting the 
continuous influx and efflux of energy between the system 
and the bath and the continuous streaming of coherence 
from the system into the bath degrees of freedom. This 
continual "production" of coherence is ultimately traced 
to the non-local nature of the quantum potential, Q. 

We present here a novel extension of the de 
Brogile/Bohm quantum theory of motion into Liouville 
space and use this to propagate quantum trajectories for 
the density matrix of a system in contact with a thermal 
bath. This approach opens a clear avenue to a number 
of new semi-classical and quantum-classical approxima- 
tions for the quantum density matrix. Furthermore, the 
formalism itself offers an interesting dynamical twist to 
interpreting decoherence and population relaxation. 
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